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Abstract 

In this investigation we revisit the question of the linear stability analysis 
of 2D steady Euler flows characterized by the presence of compact regions with 
constant vorticity embedded in a potential flow. We give a complete derivation 
of the linearized perturbation equation which, recognizing that the underlying 
equilibrium problem is of the free-boundary type, is done systematically using 
methods of the shape-differential calculus. Particular attention is given to the 
proper linearization of the contour integrals describing vortex induction. The thus 
obtained perturbation equation is validated by analytically deducing from it the 
stability analyses of the circular vortex, originally due to Kelvin (1880), and of the 
elliptic vortex, originally due to Love (1893), as special cases. We also propose and 
validate a spectrally-accurate numerical approach to the solution of the stability 
problem for vortices of general shape in which all singular integrals are evaluated 
analytically. 

Keywords: Vortex Dynamics; Linear Stability; Free-Boundary Problems; Shape 
Calculus 

1 Introduction 

In this work we are interested in the linear stability of solutions of the two-dimensional 
(2D) Euler equations which are steady in the appropriate frame of reference and feature 
compact regions with constant vorticity embedded in an otherwise potential flow. Such 
solutions arise as inviscid limits of actual viscous flows, and therefore find numerous 



applications in different areas of fluid mechanics. The first linear stability analysis of such 
a flow is attributed to Kelvin [1], see also [2, 3, 4], who identified the dispersion relation 
of neutrally-stable travelling waves with m-fold symmetry moving along the perimeter of 
a circular vortex. This problem was revisited by Baker [5] who used a contour dynamics 
approach involving integrals with singular kernels defined on the vortex boundary. The 
stability of a rotating ellipse, the so-called Kirchhoff vortex [6], was analyzed in the 
seminal study by Love [7]. The effect of an ambient strain field was investigated by Moore 
& Saffman [8], whereas Mitchell & Rossi [9] explored the relation between the linear 
instabilities and the long-term nonlinear evolution of elliptic vortices. In this context 
we also mention the analytical study of Guo et al. [10] where an integro-differential 
perturbation equation similar to ours is proposed for the case of the elliptic vortex. 
The linear stability of more complex vortex configurations, such as polygonal arrays 
of corotating vortices and translating vortex pairs, was investigated by Dritschel [11, 
12, 13], see also Dritschel & Legras [14]. A noteworthy feature of their approach is 
that they also used a continuous perturbation equation independent of a particular 
discretization employed to obtain the equilibrium solution. Problems related to the 
stability of polygonal vortex arrays were also studied by Dhanak [15]. The linear stability 
of the so-called V-states [16], corotating vortex patches and infinite periodic arrays of 
vortices was investigated in detail by Kamm [17], see as well [4]. This approach was 
based on the representation of the solution in terms of the Schwarz function, and required 
discretization and numerical differentiation in order to obtain the perturbation equation. 
A discrete form of the perturbation equation was also used by Elcrat et al. [18] in their 
investigation of the linear stability of vortices in a symmetric equilibrium with a circular 
cylinder and a free stream at infinity. There have also been a number of analytical 
investigations concerning the stability, both in the linear and nonlinear setting, of flows 
obtained as perturbations of the Rankine vortex [19, 20, 21, 22, 23, 24]. Another family 
of approaches is based on variational energy arguments going back to Kelvin. They 
were initially investigated in [11, 25, 26, 27, 28], and were more recently pursued by 
Luzzatto-Fegiz & Williamson [29, 30, 31, 32, 33, 34]. They rely on global properties of 
the excess energy vs. velocity impulse diagrams and provide partial information about 
the linear stability properties without the need to actually perform a full linear stability 
analysis. 

To the best of our knowledge, there is no complete derivation and validation of the 
stability equations for two-dimensional (2D) vortex patches available in the literature. 
The main goal of this paper is thus to fill this gap by proposing a systematic approach 
to the study of the linear stability of vortices which is quite general in the sense that it is 
based on the continuous, rather than discrete, formulation of the governing equations and 
a rigorous treatment of boundary deformations. This aspect, namely, that linearization 
of differential and integral expressions with respect to perturbations of the shape of 
the domain requires special treatment with methods other than the traditional calculus 
of variations does not seem to have received much attention in the vortex dynamics 
literature. Therefore, the proposed approach may serve as a template for studying 
more complicated problems such as the stability of vortex equilibria in three dimensions 
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(3D). Recognizing the Euler equation describing finite-area vortices as a free-boundary 
problem, we apply methods of the shape-differential calculus to derive a perturbation 
equation in the form of an integro-differential equation defined on the vortex boundary. 
We emphasize that, in contrast to [17, 18], numerical differentiation is not required at any 
stage of the derivation, but only in the evaluation of the coefficients of the perturbation 
equation. In this sense, the formulation we propose is a contour-dynamics analogue of 
the Orr-Sommerfeld equation used to study the stability of viscous parallel flows [35]. 
In this work we also demonstrate using analytical calculations how one can reproduce 
the dispersion relations obtained by Kelvin in the stability analysis of the Rankine 
vortex [1] and by Love in the stability analysis of the elliptic vortex [7] from our general 
perturbation equation. Generalization of our approach to the stability analysis of more 
complicated vortex flows, including 3D axisymmetric configurations, is deferred to future 
research. 

We are interested in developing a general approach to analyzing the linear stability of 
2D steady-state solutions of the Euler equation which is usually written in the following 
form [36] 

Aip = F(ip) in Q, (la) 
ip — ipb on dQ., (lb) 

where ip is the streamfunction, Q C M 2 is the flow domain, whereas ip b represents the 
boundary condition for the streamfunction. The function F : R -> 1 is a priori 
undetermined and any choice of F will yield a steady solution of the 2D Euler equation. 
One choice of the function F, motivated by the flow models arising from the Prandtl- 
Batchelor theory [37], is 

Fty) = -uHty -1>), (2) 

where H(-) is the Heaviside function, whereas u,i/)o EM. are, respectively, the vorticity 
inside the vortex and the value of the streamfunction at its boundary. The profile 
given in (2) corresponds to finite-area vortices with constant vorticity u embedded in 
an otherwise potential flow. This formulation can be made more general by including 
in (2) a term proportional to Dirac delta function representing a vortex sheet on the 
vortex boundary [38] , but we do not consider this in the present study. A salient feature 
of problem (l)-(2) is that the shape of the vortex region, which we will denote A, is 
a priori unknown and therefore must be determined as a part of the solution to the 
problem. System (l)-(2) is thus a free-boundary problem, and in order to emphasize 
this fact which will play a central role in our analysis, we rewrite it in the following form 
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which makes the free-boundary property more evident 



A^i = — oj 


in A, 


(3a) 


A^ 2 = 


in n\A, 


(3b) 




on dA, 


(3c) 


dn dn 


on dA, 


(3d) 


^2 = A 


on dn, 


(3e) 



where ip\ = i/j\a and ip2 = V'lnyA are the restrictions of the streamfunction ip to, re- 
spectively, the rotational and irrotational part of the flow domain and n is the unit 
vector normal to the vortex boundary ("=" means "equal to by definition"). Solu- 
tions of system (3) depend on two parameters u> and ip or, equivalently, the vortex area 
\A\ = f n H(ifj —ifj) dn and its circulation 7 = u\A\. There is evidence that by continuing 
the solutions of (3) corresponding to different flow configurations such that 7 = Const 
and \A\ — > one can obtain the corresponding point- vortex equilibria [39, 40, 41]. In 
agreement with earlier studies [17, 13], we will assume in this investigation that both 
parameters \A\ and 7 are fixed, so that the perturbations will be allowed to modify the 
shape of the vortex region A leaving its area and circulation unchanged. In other words, 
no vorticity is created by the perturbations. Derivation of the Jacobian, needed for the 
linear stability analysis of a given problem, requires differentiation of the governing equa- 
tion with respect to the state variables. In the present problem, since governing system 
(3) is of the free-boundary type, this is a rather delicate issue, because the shape of the 
vortex A can be regarded as a "state variable" . A proper way to deal with this problem 
is therefore to use the shape-differential calculus [42, 43, 44] which is a suite of math- 
ematical techniques for differentiation of various objects such as differential equations, 
integrals, etc., defined on variable domains. Application of the shape-differentiation 
methods makes it possible to derive a continuous perturbation equation for vortices 
with very general shapes in a rigorous and systematic manner which is the main con- 
tribution of this work. A key element is a proper linearization of the contour dynamics 
equations with respect to arbitrary boundary perturbations. As regards derivations of 
the different shape-differentiation identities we use in this study, the reader is referred to 
monographs [42, 43, 44], and to the review paper [45] for a survey of applications of the 
shape calculus to different problems in computational fluid dynamics. The structure of 
the paper is as follows: in the next Section we derive a general form of the perturbation 
equation, in Sections 3 and 4 we use this perturbation equation to obtain the dispersion 
relations for the circular and elliptic vortices, comments about an accurate numerical 
technique for the solution of the perturbation equation are presented in Section 5 to- 
gether with some computational results, whereas conclusions are deferred to Section 6. 
Some technical results are collected in Appendix. 



4 



2 Derivation of Perturbation Equation 



In this Section we develop a systematic approach to derivation of the perturbation equa- 
tion characterizing the stability of steady vortex configurations. For brevity of notation, 
in what follows we will use the complex number representation of vector quantities. This 
is indicated by using ordinary letters for vector quantities (e.g., n 6 C for n 6 I 2 and 
similarly for other vector quantities). Juxtaposition (i.e., Z\Z 2 for some z±, z 2 G C) in- 
dicates complex multiplication of Z\ and z 2 . However, when we need equations written 
in the vector form, as we will in referencing results from shape calculus, we will denote 
vector quantities in bold face and inner products of vectors using the usual dot notation. 
Thus, for a,beR 2 vectors and a,b e C the corresponding complex numbers, we have 
the following identity a ■ b = $R[a6], where the overbar denotes complex conjugation. 
To fix attention, we will focus on the simplest, yet generic, case of a single vortex re- 
gion A in equilibrium with an external velocity field U = u — iv which may represent 
a free stream, the influence of other vortices, or may result from transformation of the 
problem to a moving (i.e., translating or rotating) coordinate system in which a relative 
equilibrium exists. We will consider a point z = x + iy on the vortex boundary dA, see 
Figure 1, which is assumed to have counterclockwise orientation. Using techniques of 
vortex dynamics [4], the velocity of a particle located at the point z can be expressed in 
the fixed frame of reference as follows 1 

= V (t,z ) = — i dw + U(r,z), 4 

dr An J dA z — w 

where the integral on the right-hand side (RHS) has a removable singularity (w is a 
complex integration variable). Starting from relation (4), one can pursue our objective 
with one of the two alternative approaches, namely (Figure 1): 

• Lagrangian approach in which one follows the trajectory of a fluid particle on the 
vortex boundary; we note that displacement of a material particle in general also 
involves a component tangential to the vortex boundary which does not lead to 
deformations of the contour [5], 

• geometric approach in which one considers displacements of the points on the 
boundary in the direction normal to the boundary only [13]. 

Since the Lagrangian approach involves additional (tangential) degrees of freedom, the 
stability results obtained with the geometric approach represent a subset of the results 
obtained with the Lagrangian approach. While both approaches can be pursued using 
methods of shape-differentiation, in this work we will focus on the geometric approach 
as it leads to somewhat simpler calculations. 

As a first step, we transform relation (4) to a suitable moving (i.e., rotating or 
translating) coordinate system, characterized by the rigid-body velocity V (z), in which 

1 We remark that the same integral appears in [4, 47] with the opposite sign. This is due to a different 
orientation of the contour in these studies. 
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Figure 1: Schematic illustrating (solid line) the equilibrium vortex boundary dA and 
(dashed line) its perturbation. The points z e and z e represent the images of the point z 
under, respectively, the Lagrangian and "geometric" displacement. 



a relative equilibrium is attained. Denoting ((r) the position of the point z(r) in this 
moving coordinate system, the left-hand side (LHS) of (4) becomes 0- — Vo(z 
that relation (4) can rewritten as 



d( uj f z — w 
dr Air J dA z — w 



dw + U(r, z)-V (z)=V 1 {t,z) 



(5) 



In equation (5) and hereafter we use both z(r) and ((t), an d note that, since £(0) = z(0), 
either coordinate can be used on the RHS of (5) to derive the perturbation equation. 
Let n = n x + in y be the unit vector normal to the boundary dA and pointing outside of 
A. The relative equilibrium is then characterized by the condition 



dC 

d:T 



n 



u 
Air 



z — w 



OA 



W 



dw + U(t,z) - V (z) 



(6) 



which, for given vorticity oj and vortex area \A\, can interpreted as a condition on the 
shape of the vortex boundary dA. Next, we introduce a parametrization of the contour 
z = z(t, s) G dA, where r e M + is the time, such that s G [0, L] in which L denotes 
the length of the boundary dA. We will also adopt the convention that the superscript 
e, where < e < 1, will denote quantities corresponding to the perturbed boundary, 
so that ( = C\ e =o and n = n e \ e=0 are the quantities corresponding to the (relative) 
equilibrium. Then, points on the perturbed vortex boundary can be represented as 
follows 

C(T,s)=((s)+er(T,s)n(s), (7) 
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where r(t, s) represents the "shape" of the perturbation. We note that, without affect- 
ing the final result, n(s) in the last term in (7) could be replaced with its perturbed 
counterpart rf{r, s). Using equation (4) we thus deduce 



3? 



n- 



4^ 

dr 



3? 



TV 



W 



dw 



OA* 



w 



+ tt[n e (U(z*)-V (z e ))], 



(8) 



from which equilibrium condition (6) is obtained by setting e = 0. In (8) the integral 
operator is defined in terms of a contour integral. In order to be able to use shape 
calculus formulas, we need to write such integrals as integrals with respect to the arc 
length. If F is a complex-valued function and T a closed C 1 curve, we have 



<j>F(w)dw = <j> F(w)t\dw\, 



(9) 



where t is the unit tangent vector in the counterclockwise direction. 

The perturbation equation is obtained by linearizing equation (8) around the equilib- 
rium configuration characterized by (6) which is equivalent to expanding (8) in powers 
of e and retaining the first-order terms. Since equation (8) involves perturbed quantities 
defined on the perturbed vortex boundary dA e , the proper way to obtain this lineariza- 
tion is using methods of the shape-differential calculus [43]. Differentiating the LHS in 
(8) and setting e = we thus obtain 



d 
de 



3? 



n 



.d(( + em) 
~dt 





= u 


drf 


]) 






6=0 


de 



6=0 



d( d(r n) 
— + n — - — 
dt dt 



It can be shown using shape calculus [43] that 

dn e 

de 



6=0 



ds 



(10) 



(11) 



Noting also that fj!| r=0 = — (VrV) n, where Vr denotes the tangential gradient [43], 



has only a component tangential to the boundary dA and expanding 
where V lt = ft(Vi t) and §- g is 
we can rewrite (10) as follows 



dr 
dr 



dr 
dr 



dr 
Uq s i 



where V u = 3?(Vi t) and ^ is the derivative with respect to the arc length coordinate, 



d 
de 



3? 



ri 



d(( + em) 
dt 



]} 


= 3? 




6=0 



dr d( dr _ dn 
Tr t ~r + -rnn + rn — 
os dr dr dr 



dr 



(12) 



As regards the RHS in (8), we note that the dependence on the perturbation em 
is more complicated here, as the perturbation also affects the shape of the contour dA € 
on which the integral is defined. To evaluate these expressions, we need the concept 
of the "shape derivative". Given a function / : Q — > R and a perturbation vector 
field Z : Q — > IR 2 such that Z\g^ = 0, the shape derivative /' of / with respect to 
perturbation Z is defined through the identity [44] 



d 
de 



/(x + eZ(x))| e=0 = /'(x) + (V/(x)) • Z(x), V xef ,. 



(13) 
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In our problem 



/(*) 



UJ 



z 

dA z 



W 



W 



dw + U(z) - V (z) 



(14) 



and in order to apply relation (13) we need to translate the vector operations into 
complex notation. This can be done with the help of the following lemma. 

Lemma 1. For f(w, w) a complex-valued function, Z = [Zi, Z 2 ] and Z = Z\ + iZi, the 
directional derivative of f in the direction Z can be expressed as 



V/-Z 



df' df- 

fT Z+ fi= Z i 
ow ow 



(15) 



where £ 4 U£--i d 



, and JL = \{-jr + iJr) when w 

ay ' aw 2 v ox ay ' 



dw 2 V dx 

We can then write for the right hand side of (8) 



x + iy. 



d_ 
de 



{*[" , '1}L = 



d_ 
de 



= 3? 



3? 

dn e 



n 



J 8/ 



Z e — W 



de 



4tt j 8A e z c —w 
f + n(f' + r 



dw + n e (U(z e ) - V (z e )) 



e=0 



e=0 



n 



d_ 

dz 



(16) 



where we assume Z = rn in (13). 

As concerns the shape derivative /' in (16), since the reference velocity of the moving 
coordinate system Vq{z) does not depend on the perturbation, we have [Vo(^)]' = 0. For 
the integral term we need to use the following shape-calculus formula [43, page 354, 
equation (4.17)] 



g \dw\ 



a4 E 



dA 



</+(|j| + K<7)(Z.n) 



\dw\, 



(17) 



where g : Q — > K is a differentiable function (dA e C Q), g' its shape derivative and k 
is the curvature of the contour dA = T. While in the standard texts [43, 44] formula 
(17) is derived for real-valued functions g, an extension to complex-valued integrals is 
straightforward. For us, then, cf. (4) and (9) 



to z e (s) — w(q) 

g(s,q) = —^ v^t, 

4-7T z e (s) — w(q) 



s,qe [0,L]. 



(18) 



Hereafter we will adopt the convention that the subscripts z and w on different symbols, 
e.g., n, t, or r, will indicate where the corresponding quantity is evaluated. As regards 
the shape derivative g' in (17), we note that g depends on the perturbation Z through 



the independent variables only, so that by (13) we have 



dg(e,Z) 



le=0 



Vg ■ Z and therefore 



g' = 0. Finally, assembling (12) and (16), using (17), and performing the differentiations 
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required in (16) and (17), we obtain after rearranging terms the following perturbation 
equation 



dr Tr dr 
or os 



n z U'(z) + r z n z 



-— 

47T 









— w ( . 


n z (b 




JdA z 


— w\ 



dU _ dV 

dz dz 



t w \dw\ 



-n z <f> 
Jd 

f 

J 8 



+n z r 



(z — w)n w — (z — w)n v 
(z - w) 2 

n z (z — w)— n z (z — w) 

dA 



(19) 



(z — w)' 



r w t w \dw\ 



t w \dw\ 



describing the linear evolution of the local amplitude r(r, s) of the vortex boundary 
perturbation in the normal direction [cf. (7)]. In the above we have also used the identity 
t = in which follows from the counter-clockwise orientation of dA. In combination with 
an initial condition 

r(0,s)=r (s), se[0,L}, (20) 

system (19)-(20) should be interpreted as an initial-value problem in which the solu- 
tion r(t,s) is subject to the periodic boundary condition r(r, 0) = r(r,L), Vr > 0. 
We also add that some of the integrals in (19) are to be understood in the principal 
value sense. We remark that, although some terms are present in both equations, equa- 
tion (19) appears different from the perturbation equation obtained by Dritschel [13]. 
The assumption that, to the leading order in e, the perturbations do not change the 
vortex area, i.e., \A e \ = Const + 0(e 2 ) ) leads to a restriction on the admissible initial 
perturbations r . Shape-differentiating the expression for the vortex area we obtain 



f H{i) -i))dVt 
Jn 



d£l — (p rounds — (p r ds = (21) 

JA f - J JdA JdA 



which means that, in order to satisfy (at the initial instant of time) the constant-area 
condition, we need to restrict attention to initial perturbations r with the vanishing 
mean with respect to the arc length s. We add that, assuming constant vorticity u> = 
Const, condition (21) also ensures that the vortex circulation 7 remains unchanged by 
perturbations. Construction of perturbations satisfying condition (21) is discussed in 
Section 5, cf. equation (34). In the absence of the constant-area condition (21), our 
preceding analysis would need to be slightly modified, as then the shape derivative g' 
would not vanish. 

In the following two Sections we show how the stability analyses of the circular and 
elliptic vortices can be derived as special cases from general perturbation equation (19). 
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3 Stability of Circular Vortex 



In this Section we demonstrate how the well-known results concerning the stability of 
the Rankine vortex originally obtained by Kelvin [1] and reviewed also in [2, 3, 4], in 
particular, the dispersion relation for the instability waves travelling along the vortex 
boundary, can be analytically deduced directly from general perturbation equation (19). 
Even though these results are classical, we go through these calculations in some detail, 
because they will also play an important role in the stability calculations for vortices 
with arbitrary shapes. More specifically, as will be shown in Section 5 of this work, 
numerical stability computations for vortices with general shapes can be reduced to a 
"perturbation" of the corresponding analysis for the circular (Rankine) vortex. The main 
advantage of this approach is that it allows us to compute all singular (principal-value) 



integrals analytically. 

n(p) = e ip , and k = 



For the Rankine vortex with radius a we thus set z(p) 



ae 



ip 



where p = -f-s G [0, 2tv]. As regards the moving coordinate 



system, due to the rotational invariance of the problem, its velocity Vq(z) can be selected 



arbitrarily and without loss of generality we choose Vq(z) 
(19). The perturbation equation then simplifies to 



resulting in V\ t = ~tua in 



dr 



2 dp 4tt 



2tt 



.dr 
dq 



r(q) ) dq 



An 



2tt 



■ e ig + e <P) r (g) 



dq 



(22) 



where we have also written w = ae tq . In view of (21), and the periodicity of r, the first 
integral term on the RHS in (22) vanishes for area-preserving perturbations. The last 
integral in (22) corresponds to the second to last in (19) and the last integral in (19) 
vanishes (this follows from the explicit formula for principal-value integrals that we give 
below). We note that relation (22) is equivalent to the perturbation equation obtained 
with the Lagrangian approach by Baker in his analysis of the stability of the Rankine 
vortex [5]. As was done in earlier studies [1, 2, 3, 4, 5], we will adopt the "normal mode" 
approach and suppose that 



r(r,p) = p(r) e ikp + p(r) e 



j—ikp 



(23) 



where p(t) G C and k G Z. As will be demonstrated below, the normal mode approach 
is justified, because equation (22) does not lead to mode coupling. Using ansatz (23), 
the second integral in (22) becomes 

■27T ( e i q + e i«) r ( 5 ) 



Pit) 



2tt 



g i(fc+l)<2 
e ip _ e iq 
2tt e i(l-k)q 
e ip _ e iq 



dq = 
dq + e ip 
dq 



2tt 



^ikq 



e ip — e iq 



ip 



2n 



,—ikq 



e ip _ e iq 



dq 



dq 



+ 



(24) 



Each of the integrals on the RHS in (24) can be written in the form u™* iq dq for 
some M G Z. This principal-value integral can be evaluated analytically using contour 
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integration (see Appendix A) and the result is 



f 

Jo 



2* lMg f_ ne HM-l)p M>1, 



ev-e^ i 7re i(M ~ 1)p M < 1. 



Inserting ansatz (23) into perturbation equation (22) and then using (25) to reduce (24) 
we obtain an equation in which all terms are proportional either to e lks , or to e~ lks . The 
fact that there are no other terms justifies a posteriori the normal mode approach (23). 
Finally, isolating the terms proportional to e tks and e~ lks we obtain 

j t P{r) = --^k p{T) + - p{T) (26) 

and its complex conjugate copy, so that the solution is p(r) = p e^( 1_fc ) r , where po E C 
is the amplitude of the initial perturbation. From (23) we obtain the following expression 
for the evolution of the perturbation 

r(r,p) = Po A^ 1 -Q T+kp ]+p e-^ 1 -Q T+kp ] (27) 

in which the phase velocity v p = f ^i^ agrees with Kelvin's solution [1, 2, 3, 4]. We 
add that the problem analyzed by Lamb [2] is in fact formulated in a bit different way, 
because he assumed the perturbation in the form of a potential modification of the 
velocity field, so that the shape of the region A where u^O remains unchanged [2]. It 
is therefore interesting that the solutions obtained in these two problems are identical, 
despite subtle differences in the underlying assumptions. 



4 Stability of Elliptic Vortex 

In this Section we demonstrate how the well-known results concerning the linear stability 
analysis of Kirchhoff 's elliptic vortex [6] obtained initially by Love [7] can be analytically 
deduced from general perturbation equation (19). Let the vortex boundary be described 

as 

z(p) = a cos(p) + % b sin(p), p E [0, 2ir], (28) 

where a and b are, respectively, the major and minor axes of the ellipse, and p coincides 
with one of the coordinates in the elliptic coordinate system. The key result is that 
the flow becomes linearly unstable when the ellipse aspect ratio c = a/b > 3. As the 
aspect ratio increases, consecutive eigenmodes are destabilized. They are given in the 
form [7, 9, 10] 



rkip) 



dz 
dp 



cos(kp) — 



sin(kp) 



k — 3, 



(29) 
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where 



UJ 



c-1 



+ 



Ike 



UJ 



c-1 
c~+l 



Ike 



30) 

As can be verified, representation (29) satisfies the area-preservation condition (21). The 
associated eigenvalues are 



UJ 

2 



Ike 



-1 - 



c-1 

^+1 



2k 



(31) 



so that eigenmodes (29) grow proportionally to e~' tXkt . We note that only modes with 
k > 3 can become unstable, and to every such mode there corresponds a decaying mode 
with \ k = /i~ . These results are deduced from perturbation equation (19) by 

plugging r(t,p) = e~ tXkt r^(p) into this equation, where r&(p) is given in (29), and then 
performing the following steps 

1. substitute representation (28) for z = z{p) and w = w(q) in (19) and convert the 
contour integrals to definite ones with q e [0, 2n] as the integration variable, 



2. 



since p and g appear only as arguments of trigonometric functions, use the following 
identities (u G C) 



cos(kq) = 



u 



u 



sin(/cg) = 



u 



u 



(32) 



2 v " 2i 

and likewise for cos(/cp) and sin(/cp) using t> e C, to transform the resulting 
expressions to contour integrals over the unit circle |u| = 1 with du = iudp, 

3. perform partial fraction decomposition of thus obtained integrand expressions and 
evaluate the integrals using the residue theorem for the poles uq inside the unit 
circle (|w | < 1) and formula (24) for the poles u on the unit circle (|u | = 1)- 

For any given value of k, operations described in steps 1.-3. above can be performed 
symbolically using a symbolic algebra package such as Maple (the corresponding code is 
available on-line as supplementary material). Combining all resulting terms and noting 



that for the elliptic vortex Vt(p) = $l[Vi\ 
Vo(p) ■ 



ui a 2 b 2 (a sin(p) 2 +6 cos(p) 2 ) 



cf. (4), and 



-IUJ-, 



a + b ^/a^b 2 sm(p) 2 +b 4 a 2 cos(p) 2 ) ' 

z(p), cf. (5), we recover expression (31) for the eigenvalue. We add 



' (a+b) 2 

that in the special case when a = b results of Section 3 are obtained, except that in this 
limit the velocity characterizing the rotation of the coordinate system is V ot 



4 ' 



5 Accurate Numerical Solution of Stability Equa- 
tions 

In this Section we describe a spectrally-accurate numerical technique which allows us to 
solve perturbation equation (19) numerically for vortices of arbitrary smooth shape. A 
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key feature of this approach is that the singular parts of the integrals in (19) are in fact 
evaluated analytically. We remark that such an approach was also used in other stabil- 
ity analyses of 2D vortex flows, e.g., in [13]. Our method is validated by comparing the 
results obtained numerically using different resolutions for the elliptic vortex with the 
analytical results derived by Love [7], and discussed in Section 4. First, we assume that 
the vortex boundary OA admits a smooth periodic and otherwise arbitrary parametriza- 
tion z = z(p), p G [0, 2ir]. In order to simplify the notation, we rewrite perturbation 
equation (19) as 



dr(p) 



= -Mp) |w 



dr(p) 
dp 



(33) 



+ & [(A r) (p) + (£ 2 r) (p) + (£ 3 r) (p) + (£ 4 r) (p)] , p G [0, 2tt] , 



where the dependence of the perturbation r on time r was suppressed for brevity, and 
C\ r = 3? \n U'(z) + rn — ^)] , whereas £2 r , £3 r and £ 4 r represent, respectively, 
the three integral operators on the RHS of (19). Quantities appearing in equation (33) 
which are related to the contour geometry, i.e., | ^ | , t, n and k, can be evaluated using 
spectral differentiation. Equation (33) is discretized by using the following ansatz for 
the perturbation 



r(r, p) PS r (r, p) = e 



JAt 



dp 



N/2 N/2-1 

^2a k cos(kp) + ^ /3fcsin(/cp) 



k=0 



k=l 



(34) 



depending on N real-valued coefficients a k and (3k and where A G C is the growth rate 
we seek to determine. The factor | ^ | in (34) ensures that area-preservation condition 
(21) is satisfied by all terms in the series, except for the constant term a cos (0p) which 
must be introduced in order to ensure well-posedness of the collocation problem (it will 
be, however, subtracted via projection from the resulting eigenvectors). After using the 
ansatz (34), the independent variable p is discretized using N equispaced collocation 



points {pj}j =1 (N is an even number). Denoting y = [a , . . . , ajv/2, 
W N , this discretization leads to the following eigenvalue problem 



Ay = -^A- 1 (B + C)y, 



where 



— fiipj), j,l = l,...,N, 



cos((/ - l)p), 



in which n(p) 



sin 



is the (invertible) collocation matrix, 



N 



N 

I = !,...,- + !, 
N 

, l = - + 2,...,N 



^N/2-l] T £ 



(35) 
(36) 
(37) 



^ 1—1 , 

i B h = -Vitipj) j-(Pj)\ 5^A jm D mi , j,l = l,...,N 

m=l 



(38) 
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corresponds to the first term on the RHS in (33) with D denoting the spectral differen- 
tiation matrix associated with ansatz (34), and 



[Cfc, = HKArOfo) + (An)fe) + 0C 3 r,)(p,-) + (An)fe 



j,/ = l,...,iV (39) 



in which the first term on the RHS is obtained via straightforward evaluation whereas 
the calculation of the integral expressions C^Ti, C$ri, and £4^ is described below. 

All integrals with bounded integrands are evaluated using the trapezoidal rule, which 
on a periodic domain is a spectrally-accurate quadrature. Since p and q are now the 
independent variables, in this Section these characters will be used as subscripts to 
indicate where a given variable is evaluated. The integrand expressions in V u , cf. (4), 

and in C 2 r 4 1 J - K r g ) i, ||| dq have removable singularities 

and are therefore bounded [4]. As regards the integral operators £3 and £4 which 
are defined in the principal- value sense, we use a standard approach (see, e.g., [48]) 
in which they are split into principal-value integrals corresponding to a circle, which 
can be evaluated analytically using formula (25), and regular "corrections" which can 
be approximated numerically with the spectral accuracy using the trapezoidal rule. In 
order to be able to use (25), the basis functions in ansatz (34) are expressed in terms of 
the complex exponentials, so that 



£3 
£3 



dz -J 
dq 

dz - ! 



dq 



cos{k q) 
sin(/c q) 



£3 
£3 



dz 



-1 



Jkq 



dq 

dz - 



dq 



ikq 



+ £3 



dz 



-1 



dq 

dz -1 



dq 



-ikq 



-ikq 



and analogously for £4 r. Then, we have 



£,, 



dz 



dq 



,;irnq 



(P) 



u n 



47T 



dz 
dq 



(?) 




2w J^± n p- n v dw 



z p — w q dq 



dq 



,jimp 



Q(p, q)dq + i 



2tt e iq 



e ip _ e tq 



-Q(P, <l)dq> , 



where Q(p, q) 



(40) 

n p — n p . Expanding this expression in a Fourier series with respect 



to q, i.e., Q{p,q) = J2k=-N/2 Qk{p) etkq i the second integral on the RHS in (40) can be 
evaluated as follows with the help of formula (25) 



,2* et q W a , 



k=-N/2 
-1 



= VK 



"2tt e i(k+l)q 
e ip _ e iq 



N/2 



dq 



Qk(p)e ikq -J2Q^P) elkq 

k=-N/2 



k=0 



(41) 
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As regards the remaining integral operator, we have 



dz 
dq 



-i 



D vmq 



n q dw 




dz - 1 

— (q) S{p,q)dq 



e «p _ e iq 



dz - 1 

— (q) S(p,q)dq}, 



(42) 

where S^p, q) = Ti q — n q . Expanding this expression multiplied | j^(q) | in a Fourier 

series with respect to q, i.e., |^(?)| 1 S(p,q) = J2k=-N/2 (\%\ ^G )^ e%k9 i the second 
integral on the RHS in (42) can be evaluated as follows with the help of formula (25) 



1-2-K 

Jo 



tmq 



dz 
dq 



(?) 



2vr e i(k+m+l)q 




111 



The first integrals on the RHS in (40) and (42) are regular and can be evaluated nu- 
merically (they also vanish in the circular case). Expressions £3 ^|^| 1 e~ tmq ^j and 

£4 ^|^| 1 e~ vmq ^j are computed in an analogous manner. After collecting expressions 

(36)-(43), the algebraic eigenvalue problem (35) can be solved using standard tools. 

We close this Section with a numerical validation of the proposed method. As the test 
problem we used the case of the elliptic vortex, also studied analytically in Section 4. We 
compute the eigenvalues by solving problem (35) using different numerical resolutions 
iV and compare the results against the closed-form expression (31) obtained by Love 
[7], see also [10, 9]. More specifically, we analyze the imaginary parts 3(A) of the 
eigenvalues responsible for the instability. In Figure 2 we show the relative errors between 
the eigenvalues Q(X N ) computed numerically with resolution N and ^s(Xk) obtained 
analytically by Love, cf. (31), for two different aspects ratios of the elliptic vortex (for 
the aspect ratio equal to 8 there are in fact 4 distinct unstable eigenvalues). In Figure 
2 we note that this error decays exponentially fast dropping to the machine precision 
level already for modest resolutions N, thereby confirming the spectral accuracy of the 
proposed method. We add that, as expected, the resolutions required to achieve a given 
accuracy increase with the aspect ratio of the vortex. 
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[ , , , , jg , , , , , , i 

50 100 150 200 250 300 

N 

Figure 2: Relative errors in the imaginary parts of the numerically computed eigenvalues 
X N with respect to the analytical expressions for given in equation (31), cf. [7], as 
a function of the resolution iV for vortices with (open symbols) aspect ratio a/b = 4 
and (solid symbols) aspect ratio a/b = 8. For the aspect ratio a/b = 4 the eigenmode 
wavenumber corresponding to the eigenvalue is k = 3, whereas for the aspect ratio 
a/b = 8 the wavenumbers are k = 3,4,5,6, cf. (34), with larger values of k resulting in 
larger errors. 

6 Conclusions 

While perturbation equations similar to (19) have already been used for the linear sta- 
bility analysis of 2D vortex patches [11, 14, 13], to the best of our knowledge, the present 
study offers a first complete derivation of this approach. In particular, it relies on meth- 
ods of the shape calculus which are a general and mathematically consistent way of 
dealing with the free-boundary aspect of the problem. We add that, in the context 
of vortex dynamics, such techniques have already been used to study continuation of 
families of solutions [41] and vortex control problems [49]. The proposed numerical ap- 
proach is demonstrated to be spectrally-accurate reducing all singular integrals to closed 
form (24). The two validation tests presented offer "the proof of the concept" for this 
approach. 

Generalization of our method to problems involving several vortices is straightfor- 
ward, and requires that an equation of the form (19) be written for each individual 
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vortex with the interaction between the vortices captured by the field U (z) (which van- 
ishes identically in the single- vortex examples considered in Sections 3 and 4). This 
description will be simplified by symmetries of the vortex configuration. One aspect 
of this problem which does not seem to have received much attention in the literature 
is the possibility for subcritical disturbance amplification due to non-normality of the 
underlying stability operator. While these questions have already been investigated for 
viscous vortices with unbounded vorticity support (e.g., [51, 52]), to the best of our 
knowledge, there are no results concerning free-boundary problems of the type (l)-(2). 
Another interesting and far less researched application area for our approach is the sta- 
bility analysis of 3D axisymmetric vortex flows, with and without swirl, with compact 
vorticity support [50] . All these problems are left to future research. 

As regards limitations of the proposed approach, we remark that shape calculus in 
its standard formulation [46] is applicable only to problems posed on smooth manifolds, 
hence our method would need to be modified, so that it can be applied to contours with 
singularities, such as, e.g., corners (flows with such features typically arise as terminal 
members of families of steady solutions [41]). Some relevant ideas are already mentioned 
in [53]. Likewise, analysis of stability with respect to perturbations leading to topology 
changes requires the use of different methods. 
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A Evaluation of the Integral f^ n e f s ^ iq dq 

The calculations presented in Section 3, 4 and 5 required the values of the integrals 

/•27T AM p2w iM 

1= / -dq = e~ ip / 7 ,dq (44) 

for integer M. If M > 0, let us consider 

e~ ip / — dq 

in which C is now the contour connecting the points 0, 2%, 2% + Yi, Yi, and with 
a small semi-circular indentation into the upper half-plane above p, where p G (0,2%). 
This contour integral vanishes by Cauchy's Theorem. The integrals over the left and 
right lateral sides of C cancel by periodicity. As regards the top segment of contour C , 
writing q = x + iY , we have 

e iMq = e iMx e -MY ^ jiq-p) = e i(x~p) & -Y ? 
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so for Y — > oo the integral over that segment goes to zero. There remains to take 
the limit of the integral over the bottom segment of the contour C as the radius of 
the semi-circular indentation vanishes. This is the original principal-value integral (44) 
augmented by the contribution from the indentation, which can be evaluated using the 
Cauchy integral formula. To this end, we write the integrand as 

e iMq (q - p) 1 

1 _ e i(q-p) q -p 

and find, using L'Hopital's rule, 

lim e^fa-rt =feitM _ 1)r 
q^P 1 - e l (v-P> 

Thus, the limit of the integral over the indentation is this value multiplied by —iri, and, 
finally, we obtain 

/ = 7 ri(ie i < Af - 1 > p ) = -ne* M -V p . (45) 

When M < 1, instead of contour C, we use its reflection into the lower half-plane. Then, 
since Y is negative, the integral over the bottom segment goes to zero as Y — > — oo and, 
since the small semi circle is positively oriented, 

/ = - n i(i e W-VP) = 7re <(^-i)f. (46) 
Formulas (45) and (46) are equivalent to (25). 
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